# Decision analytic model {#model}

```{r setup2, include=FALSE, echo = FALSE, message = FALSE, warning=FALSE}
knitr::opts_chunk$set(echo=TRUE)
knitr::opts_chunk$set(message = FALSE)
rm(list = ls())  # delete everything that is in R's memory
options(scipen=1000) # removal of scientific notation

path <- "C:/Users/PouwelsXGLV/Documents/Maastricht/Eribulin CEA/Model final_Gamma/" #set location of model, to adapt!

setwd(paste(path,"report", sep = "")) # set working directory to File Source

#load package
library(rms)
library(survival)
library(muhaz)
library(survminer)
library(ggplot2)
library(data.table)
library(printr)
library(dplyr)
library(flexsurv)
library(tidyverse)
library(scales)
library(matrixStats)
library(reshape2)

# load functions
source(paste(path,"R/01_fun_analyses.R", sep = "")) 
# load model inputs
load(paste(path,"data/01_model_inputs.R", sep = ""))
```

This section presents the cost-effectiveness models used for the analyses. The current cost effectiveness analysis was composed of two main functions. The `CEModel_det` function which estimates the deterministic results and the `CEModel_prob` which estimates the probabilistic results. Both functions are defined in the *01_fun_analyses.R*  (`R` folder). Each function only requires the user to define the list of parameters to use in order to obtain deterministic and probabilistic results. In the case of the deterministic model, the `params_det` list can be used (saved in the *01_model_inputs.R*  file), for the probabilistic model, probabilistic parameters have to be estimated first.

## Deterministic model
In the deterministic model, the estimated survival curves for PFS and OS are directly used to determine the state occupancy of the PF, PD, and death health state. The TTD health state is used to estimate the chemotherapy costs accrued in the PF health state.  
In our analysis, the `CEModel_det` function is informed by the `params` argument. This argument determines the model inputs to use during the deterministic analysis. To obtain the deterministic results presented in the manuscript, we have created the `params_det` list of parameters, which is also saved in the *01_model_inputs.R* file (`data` folder).  
The deterministic model is divided in three parts: the initialisation phase, the processing phase, and the output calculation phase. This model structure has been inspired by the Decision Analysis in R for Technologies in Health publication [@Alarid-Escudero2019].  
During the initialisation phase, multiple matrices and vectors that are used in the processing and ouput calculation phases are declared. In the processing phase, the cohort traces (or state occupancy) for the eribulin and non-eribulin arms are determined (respectively `m.TR` and `m.COMP` objects). Finally, aggregated and disaggregated outcomes are calculated in the output phase. The result of the deterministic model is a list (`l.det.res`) containing the aggregated (`l.det.res[[1]]`) and disaggreggated results(`l.det.res[[2]]`), as well as the state occupancy matrices of the eribulin (`l.det.res$Trace.trt`) and non-eribulin group (`l.det.res$Trace.comp`). In the *01_model_inputs.R* file, the `CEModel_det` function is extensively anotated (`R` folder).

```{r, echo = TRUE}
print.function(CEModel_det)
```


```{r, echo = T}
l.det.res <- CEModel_det(params = params_det) # run the function
```

To inspect the cohort traces of the deterministic model, one can use the following code. 
```{r, echo = T}
head(l.det.res$Trace.trt[,1:3]) # start of the eribulin cohort trace
# the fourth column is the TTD column, which is only used for the calculation of the systemic treatment costs.
```

By using the code below, one can inspect the cohort traces of both treatment arms.
```{r fig-traces, echo=T, fig.cap = "Cohort traces"}
## create time variable
Time <- c(0:abs(n.tw))*n.days.week/(365.25/12)
l.det.res$Trace.trt <- cbind(l.det.res$Trace.trt, Time)
l.det.res$Trace.comp <- cbind(l.det.res$Trace.comp, Time)

plot(x = l.det.res$Trace.trt[,"Time"], 
     y = l.det.res$Trace.trt[,"PFS"],
     type = "l",
     lwd=2, col="blue", xlab = "Time (in months)", 
     ylab = "Propotion of patients in each state",
     ylim = c(0,1),
     cex.lab=0.75, main = "Eribulin") 
lines(x = l.det.res$Trace.trt[,"Time"], 
      y = l.det.res$Trace.trt[,"PD"],
      type = "l",
      lwd=2, col="red")
lines(x = l.det.res$Trace.trt[,"Time"], 
      y = l.det.res$Trace.trt[,"Dead"],
      type = "l",
      lwd=2, col="black")
legend("topright", legend = c("PFS", "PD", "Dead"), col=c("blue", "red", "black"),  
       lty=c(1,1,1), cex=0.8)

plot(x = l.det.res$Trace.comp[,"Time"], 
     y = l.det.res$Trace.comp[,"PFS"],
     type = "l",
     lwd=2, col="blue", xlab = "Time (in months)", 
     ylab = "Propotion of patients in each state",
     ylim = c(0,1),
     cex.lab=0.75,  main = "Non-eribulin") 
lines(x = l.det.res$Trace.comp[,"Time"], 
      y = l.det.res$Trace.comp[,"PD"],
      type = "l",
      lwd=2, col="red")
lines(x = l.det.res$Trace.comp[,"Time"], 
      y = l.det.res$Trace.comp[,"Dead"],
      type = "l",
      lwd=2, col="black")
legend("topright", legend = c("PFS", "PD", "Dead"), col=c("blue", "red", "black"),  
       lty=c(1,1,1), cex=0.8)
```

## Probabilistic model
The probabilistic model has been built based on the deterministic model, and consists of a loop over the probabilistic values of the effectiveness, utility, and resource use and costs parameters. The probabilistic parameters and results will be calculated in the Section \@ref(analysis) and will be  combined in a list that will be used to produce probabilistic results by the `CEModel_prob` function.  
```{r, echo = TRUE}
print.function(CEModel_prob)
```
